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In this work we consider a system of spinless fermions with nearest and next-to-nearest neighbor repul- 
sive Hubbard interactions on a honeycomb lattice, and propose and analyze a realistic scheme for analog 
quantum simulation of this model with cold atoms in a two-dimensional hexagonal optical lattice. To this 
end, we first derive the zero-temperature phase diagram of the interacting model within a mean-field the- 
ory treatment. We show that besides a semi-metallic and a charge-density-wave ordered phase, the system 
exhibits a quantum anomalous Hall phase, which is generated dynamically, i.e. purely as a result of the re- 
pulsive fermionic interactions and in the absence of any external gauge fields. We establish the topological 
nature of this dynamically created Mott insulating phase by the numerical calculation of a Chern number. 
Based on the knowledge of the mean-field phase diagram, we then discuss in detail how the interacting 
Hamiltonian can be engineered effectively by state-of-the-art experimental techniques for laser-dressing 
of cold fermionic ground-state atoms with electronically excited Rydberg states that exhibit strong dipolar 
interactions. 
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I. INTRODUCTION 

Quantum systems and materials with topological prop- 
erties and quantum order are currently the focus of in- 
tense theoretical and experimental research: On the one 
hand, they introduce a new paradigm in condensed mat- 
ter physics challenging the standard model of the Landau- 
Ginzburg theory of phase transitions based on the standard 
symmetry breaking mechanism. New topological quan- 
tum systems exhibit non-local order parameters and re- 
quire a different framework than the Landau-Ginzburg the- 
ory for their understanding and classification 1 . On the other 
hand, topological quantum systems offer new fascinating 
potential applications like topological quantum informa- 
tion processing^]. 

These developments have not gone unnoticed by an- 
other recently emerging field: quantum simulations with 
cold atoms in optical lattices. Originally conceived as plat- 
forms to unveil the complicated physics of strongly corre- 
lated quantum systems in condensed matter 4 5 , its range of 
applications has rapidly increased 6 10 . The topic has been 
boosted by enormous progress in the field of cold quantum 
gases and quantum simulation, where these systems offer 
unique possibilities to realize and study in a controllable 
way the physics of some of these models over a wide range 
of parameters. 

A particular fascinating class of new topological quan- 
tum systems are topological insulators (TIs) that provide 
a new paradigm of how topological properties can be 
even realized in fermionic free systems under appropri- 
ate conditional"^! New time reversal TIs were theoreti- 
cally proposed in two-dimensional hexagonal lattice^^ 
and soon extended to three-dimensional arrays^^. Re- 
markably one -dimensional versions of TIs can be also 
explicitly realized by means of lattice gauge theory in- 
spired models 24 " 26 . In a series of ground-breaking experi- 



ments, the existence of this new quantum phases of mat- 
ter has been unambiguously confirmed thereby establish- 
ing the topic of TIs as a solid and consolidated research 
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There are two main avenues to realize TI phases: 
i/ The original approach to simulate TI systems with cold 
atoms in optical lattices stems from the idea of synthesiz- 
ing gauge magnetic fields in an effective way 33 37 , which 
has recently been demonstrated experimentally 3 ^ 9 -'. The 
reasoning behind these quantum engineering ideas is that 
background fixed classical gauge fields play a fundamen- 
tal role in the theoretical construction of topological in- 
sulating systems. As an example of the broad range of 
possibilities of this new line of quantum simulation com- 
pared to the study of standard condensed matter systems 
in nature, it is worth mentioning the possibility of simulat- 
ing non-Abelian gauge fields in several different cases^ ! 4 ]!, 
in particular an anomalous relativistic non-Abelian quan- 
tum Hall effect 43 . These synthetic gauge field implemen- 
tations have paved the way to the quantum simulation of 
TIs with optical lattices, ranging from proposals for one- 
dimensional systemP^ElllIS two-dimensional l attic esP^^ 
even up to three-dimensional implementations 57 58 involv- 
ing realizations of lattice gauge theory effects. 

ii/ Alternatively, there exists the less explored possibil- 
ity that fermionic interactions lead to the dynamical cre- 
ation of a TI phase: Understanding the role of fermionic 
interactions in TIs is attracting a lot of attention since the 
classification of non-interacting TIs and superconductors 
has been achieved in the form of the celebrated "Periodic 
Table' 59 60 . This classification is exhaustive for quadratic 
fermionic Hamiltonians and depends on the dimensional- 
ity of the lattice fermions and the discrete anti-unitary sym- 
metries which are broken or not by the dynamics, specif- 
ically time-reversal and particle-hole symmetry as well as 
chirality. Interestingly enough, with the pioneering work 
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by Raghu et al. 61 and subsequent 

studiepMl 

it has be- 
come clear that interactions in purely fermionic Hamilto- 
nians may not only change the topological class of a TI into 
a different (possibly trivial) one. Instead, the interactions 
may actually become a resource to produce TI phases. The 
key point of this new mechanism is that there is no need to 
resort to background-fixed gauge fields coupled to fermion 
degrees of freedom, but that it is the very presence of many- 
body interaction effects among fermions hopping in partic- 
ular lattice geometries that produces a new phase called a 
Topological Mott Insulator (TMI). 

In this work we are interested in the latter interacting sce- 
nario iil as a way to produce TI phases whose non-trivial 
topological properties arise purely as a consequence of the 
fermionic interactions, as opposed to the scenario i/ of 
gauge-field induced TIs. We will consider a model of inter- 
acting spin-less fermions on a hexagonal lattice, which we 
specify below, where this mechanism may result in the for- 
mation of a TMI phase, generated dynamically by repulsive 
Hubbard-type interactions and in the absence of any (syn- 
thetic) external gauge fields. Besides a theoretical analysis 
of the model, we will propose how this new physics could 
be observed in an experiment by implementing the inter- 
acting fermion model in a quantum simulator using cold 
Rydberg atoms in optical lattices. On the one hand, such 
quantum simulation is important for benchmarking theo- 
retical predictions with experimental observations, in par- 
ticular as the problem of interacting lattice fermions in two 
or higher spatial dimensions is notoriously hard to solve. On 
the other hand, as we will show below, the flexibility to tune 
the parameters of the model of interest over a wide range 
of parameters in a quantum simulation offers the possibil- 
ity to explore in the laboratory the physics in regimes, which 
are not naturally realized or not readily accessible in known 
materials. 

In the context of quantum simulation, we will see that as 
an alternative to engineer complex next-to-nearest neigh- 
bor hopping dynamics as required in the Haldane model 14 , 
which supports a quantum anomalous Hall (QAH) phase 
that has so far not been observed in any material, intro- 
ducing long-range interactions between the fermions might 
become an ally, as long as these turn out to favor the for- 
mation of a stable QAH phase and there is no need to im- 
plement synthetic gauge fields. We anticipate at this stage 
that the fact that both the Haldane model and the interact- 
ing model considered by Raghu etal. 61 and also by us in this 
work are defined on the same type of hexagonal lattice is 
relevant for the following reason: the next-to-nearest neigh- 
bor fermionic interactions, when treated within a mean- 
field theory (MFT) approach, give rise to an order parame- 
ter that resembles the single-particle next-to-nearest neigh- 
bor hopping in the Haldane model. This connection is the 
physical key mechanism underlying the dynamically cre- 
ated, i.e. interaction-induced TMI phase described above. 

Finally we remark that the hexagonal lattice geometry 
considered in this work and our results are also related to 
the very rich physics in graphene 67 69 , where gauge fields 
can be induced by means of disclinations or strain fields 70 . 



A. Structure of the paper 

In this work we 

1. identify the relevant order parameters and determine 
the complete zero-temperature phase diagram within 
a MFT treatment, and compare our findings to re- 
lated, previous work on this model 61 l 62 i 

2. show that, within the considered MFT approach, 
there is no region of the phase diagram with coex- 
istence of the encountered topologically non-trivial 
Mott insulating phase and any of the other phases 
with local order parameters, 

3. underpin the character of the topologically trivial and 
non-trivial phases by an explicit numerical calcula- 
tion of Chern numbers 71 77 , which we compare with 
the behavior of local order parameters, 

4. propose and analyze a realistic scheme for an analog 
quantum simulation of the model with cold fermionic 
Rydberg atoms in an optical lattice. 

This paper is organized as follows: in Sect.[IT]we introduce 
the model Hamiltonian describing spinless fermions on a 
hexagonal lattice subject to Hubbard interactions. Then, 
we go on to review the basic properties of the free hopping 
Hamiltonian part, in particular the momentum represen- 
tation of quantum states in the Brillouin zone that consti- 
tute the basis for the rest of the analysis. 
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In Sect. 

perform a MFT analysis for dealing with the Hubbard in- 
teraction terms, which allows us to identify and character- 
ize the relevant phases in the various parameter regimes. In 



Sect. IV we show how the model Hamiltonian of interacting 



fermions can be engineered and arises as an effective de- 
scription of the dynamics of fermionic ground state atoms 
laser-coupled to Rydberg states, which exhibit strong and 
long-range interactions. Finally, Sect|V| is devoted to con- 
clusions and an outlook. 

Two appendices contain the detailed and basic explana- 
tions of the constructions used throughout the text. Specif- 
ically, Appendix|X|presents the details on the calculation of 
the Chern numbers characterizing the topological Mott in- 
sulating phase deduced from the MFT analysis. AppendixjB] 
summarizes details on various parametrizations of the Bril- 
louin zone and domain of integration. 



II. THE MODEL 

Let us consider spinless fermions at sites j e A of a hexag- 
onal lattice in two dimensions. The associated fermionic 
creation c, and annihilation operators cj satisfy the canon- 
ical anti-commutation relations {c ; -,ct} = <5,- ,-. The model 

J ' J 

Hamiltonian H we are interested in is an extension of the 
free fermion hopping H t Hamiltonian between nearest- 
neighbor sites of the lattice by including Hubbard-type 
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FIG. 1. Hexagonal lattice A with sites j decomposed in two-site 
basis cells /' = {m,n) connected by basis vectors ai and a2- A 
fermion at the 0-type site (m, n) possesses three nearest-neighbor 
sites (shaded black filled circles) and six next-to-nearest neighbor 
sites (shaded white open circles). 



nearest-neighbor H^n and next-to-nearest neighbor Hnnn 
interactions: 



H:=H t + _Hnn + f^NNN 
— -tj^tcjcj + h. 



c) + V\ E niUj + 



V 2 Y, ntnj, (1) 



where denotes the summation over pairs of nearest- 
neighbor, and over the next-to-nearest neighbor 
sites, and n,- = etc,- denotes the number operator on lattice 
site i. The real-valued hopping amplitude is t, while the in- 
teracting couplings are repulsive V\, V 2 > 0. 

This model has appeared in different contexts, such as the 
study of antiferromagnetic phases using quantum Monte 
Carlo simulations 78 , and as a candidate for a quantum spin 
liquid 79 80 . Here, our focus lies on a quantum simulation of 
the model with Rydberg atoms in optical lattices. To con- 
nect with a physical implementation in a cold atoms setup, 
we will first study the model theoretically which allows us to 
establish its phase diagram and identify the relevant order 
parameters and observables to be measured in an experi- 
ment. 

Before dealing with the full Hamiltonian [TJ we start first 
by reviewing the properties of the free part of the Hamilto- 
nian H t , in order to establish the basic notations and oper- 
ators that we shall use below when treating the interaction 



terms in Sect III The bipartite nature of the hexagonal lat- 
tice makes it convenient to rewrite the Hamiltonian H t in 
terms of two-site basis cells^ (see Fig. [TJ , labeled by an in- 
dex pair [m, n), each of which contains one site of the <p- 
sublattice (open white circles) and one of the i//-sublattice 
(filled black circles) : 



E «' 



y/m 



n c <pmn + c ymn+l Cl l> mn + c ym-ln Cl P m '^ + 



Correct counting of terms, thus avoiding possible double 
counting errors, is essential, in particular to correctly repro- 
duce the quantitative details of the mean-field phase dia- 
gram of {TJ , such as the values of the critical couplings where 
phase transitions take place. We will comment on some 
controversy on this issue in the recent literatur e 1 61 ! 62 ^ more 
detail below. At this point we would like to note that the sum 
over n, m of the hopping processes of a fermion at a site cp- 
type lattice site {m, n) to one of its three neighboring i/'-type 
sites already covers all possible hopping processes along all 
links connecting neighboring sites of the hexagonal lattice. 

In addition to the tight-binding term, we consider a stag- 
gering potential, describing a chemical potential difference 
for fermions residing at (p- and t//-sites, which reads 



c il>mn c ymn) ■ 



(3) 



To determine the energy spectrum of the tight-binding 
Hamiltonian in combination with the staggering potential, 
we introduce Fourier-transformed fermionic operators, 



1 



E exp(z'k (k), 

N A keBZ 



(4) 



with a e {(p,if/}, to rewrite the Hamiltonian in momentum 
space as 

H t + H p ^Y^^)[_f Ak "^)*(k). (5) 

Here, T^fk) = (c^fk), cJ,(k)J, and r mn = na\ + m&2 are the 
real-space vectors of the lattice, with 



ai = (3/2, v/3/2), a 2 = (-3/2, v / 3/2) 



(6) 



denoting the basis vectors (for a unit length lattice spac- 
ing), which span the whole lattice in the two-site basis 
parametrization (see Fig. [TJ . The total number of two-site 
basis cells is Na, which is half of the total number of lattice 
sites, A/a - NI2. The function 



- 1 + exp(zkai) + exp(- ika2) 



(7) 



contains the information about the structure and symmetry 
of the hexagonal lattice. The reciprocal lattice vectors, sym- 
metry properties and various equivalent parametrizations 
of the Brillouin zone, either in a discretized or in a continu- 
ous form, are provided in AppendixlB] 

Diagonalization of Hamiltonian j5j readily yields the two- 
band energy spectrum 



E ± Qt) = +dp 2 + t 2 \A k \ 2 



(8) 



(2) 



= +Y^ 2 + f 2 (3 + 2cos(v / 3fc y )+4cos(v / 3fc y /2)cos(3fc.J2)j 

Whereas a non-vanishing staggering potential /5 opens an 
energy gap between the two bands, for p = the gap 
closes and the spectrum exhibits the well-known band 
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structur^HEH with a linear dispersion relation in the vicin- with 
ity of the Dirac cones (see Fig.[9]in AppendixjBj. Our moti- 
vation to consider the staggering potential Hp at this stage 
is because interacting terms Hnn and Hnnn, when treated 
within a MFT approach in the following section, will pro- 
duce in momentum space a very similar structure. In fact, 
we anticipate that - similar to the effect of an arbitrarily 
weak staggering potential - repulsive nearest- and next-to 
nearest-neighbor interactions above a critical strength can 
induce openings of an energy gap, associated to first- and 
second-order quantum phase transitions. 



III. MEAN-FIELD ANALYSIS AND PHASE DIAGRAM 

A. Mean-field treatment for nearest-neighbour interactions 

As a first step, and for clarity of the presentation, let us 
for the moment consider the mean-field treatment of the 
Hamiltonian {!] only including nearest-neighbor interac- 
tions, 



Hi :- H t + Hnn- 



(9) 



This analysis presents some of the basic ingredients that will 
later be used in the study of the complete Hamiltonian, and 
introduces the charge-density-wave (CDW) order parame- 
ter as one of the basic order parameters of the system. 

We perform the mean-field approximation of the inter- 
acting part, which is quartic in fermionic operators, accord- 
ing to the Wick theorerrP^ 

nirij ^{c\ci)c\cj + {c\cj)c\ci - (c\ct)(c\cj) 

- {(c\cj)c]ci + {c]ci)c\cj - {c\cj){c]c t )] . (10) 

By introducing the self-consistently defined mean values 
Hi :- (etc,-) and := (etc/) the Hamiltonian H\ in mean- 
field approximation takes the form 

HlMF = " £ [(f+Vlftf)c}c* + /l.C.-Vi|{, / | 2 ] 

(11) 



(hi) 

+ Vi £ ("i n j + Hj rii - HiHj) . 

One sees that the presence of H^n leads in mean-field ap- 
proximation, besides the repulsive density-density-type in- 
teraction terms, to a renormalization of the bare hopping 
amplitude t, 



(12) 



Similar to the rewriting of H t above, we now reformulate 
H\mv in the two-site basis, and in addition assume that the 
expectation values of the fermion density on <p- and on y/- 
sites of the lattice, and H^, are translationally invariant, 
and that £ := (ctcy) is real-valued and rotationally invari- 
ant. We thus obtain the mean-field Hamiltonian in momen- 
tum space 

Himf = £* + (k)^iMF (k) * (k) - 3iV A Vy (n^ - { 2 ) (13) 



Nearest-neighbor repulsive interactions are known to break 
the symmetry of cf>- and i//-lattice sites and to favor the for- 
mation of a CDW 61 . We define the CDW order parameter 



p := - {n^-ny,) 



(15) 



which quantifies the fermion density imbalance at <p and y/- 
sites of the two sub-lattices, for a total density of 



1 _ 

n:= -(fiQ + ny) 



(16) 



in the system. Diagonalization of (14) yields the energy- 
momentum dispersion relation 



£±(k) = 3Vin± J(3Vip) 2 + (f+Vif) 2 |A k | 2 



(17) 



At zero temperature, the lower energy band is completely 
filled, and the free energy (per two-site basis cell) is 

FIN A = ?>V l {p 2 + ?)-?>V l n 2 +\ f d 2 kE_(k). (18) 

Lr JbZ 

Here, we have converted the discrete sum over momenta 
into an integral over the first Brillouin zone, with L 2 denot- 
ing the integration area in momentum space (see Appendix 
|B]for details). 

In which quantum phase the system is, as a function of 
the coupling parameters V\ and t, is determined by the 
global minimum of the free energy functional. Evaluation 
of the partial derivatives dF/d£ = and dFldp = yields the 
set of coupled self- consistency equations for the CDW order 



parameter p and the renormalization 12 of the nearest- 
neighbor hopping amplitude 



P = 



(r+<Vi)|A k | 2 



6L 2 Jbz ^/{3V 1 p) 2 + {t + ^V 1 ) 2 \A k \ 2 ' 
2L2 Jbz VoVipjZ-Kr + fVi) 2 !^! 2 



(19) 
(20) 



In the discussion of the phase diagram below, these equa- 
tions will serve us to study the renormalization of the hop- 
ping amplitude and to determine the critical interaction 
strength V\ c (for vanishing next-to-nearest neighbor inter- 
actions V2) at which a phase transition from a SM region to 
CDW phase takes place. 



B. Mean-field analysis of the next-to-nearest neighbor 
interactions 

In this subsection we analyze the effect of the Jf2-term in 
a mean-field treatment, defining the auxiliary Hamiltonian 



H2 :- H t + -Hnnn- 



(21) 
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As for the nearest-neighbor interactions, we decouple the 
quartic interaction terms according to the mean-field ap- 
proximation (TO) . 

The diagonal part (first line of (To) ) leads to an effec- 
tive next-to-nearest neighbor density- density interactions 
of fermions residing at sites of the same sublattice a e {<p, y/}. 
Inspection of one representative two-site basis cell {m, n) 
shows that an a-site possesses six next-to-nearest-neighbor 
a-sites. Each of these six next-to-nearest neighbor inter- 
action terms produces a mean-field contribution {2n a n a - 
~n 2 a ), which, however, is to be counted only with weight 1/2 
to avoid double-counting, when the sum over all two-site 
basis cells of the lattice is carried out. In consequence, after 
performing the transformation to momentum space the re- 
sulting density-type contribution to the Hamiltonian reads 



k a 



(22) 



The off-diagonal terms (second line of 10 ) play a 
fundamentally different role: They correspond to 
emerging next-to-nearest neighbor hopping processes, 



Xa,m,n,m' ,n' c am 'n' Camn 

mean values Xa,m,n,m',n 



with the self- consistently defined 



: = (CamnCam'n')- A mean-field 



ansatz for the %a mean values on the two sublattices, which 
(i) is translationally invariant, i.e. under displacements ac- 
cording integer multiples of the real-space lattice vectors ai 
and a 2 , (ii) which still possesses a rotational symmetry un- 
der rotations about an angle of 27i73, and (iii) is compatible 
with complex- valued ^-fields, is given by the identification 
(see Fig. [2} 

Xa = Xa,m+\,n,m,n — Xa,m,n+l,m,n — %a,m—\,n—\,m,nt (23) 
Xa = Xa,m+l,n+l,m,n — Xa,m-\,n,m,n — Xa,m,n-l,m,n- (24) 

It will be convenient to parametrize the two mean fields in 
terms of their absolute value and phase, 



Xa — \Xa\ e 



a e {(p, if/}, 



(25) 



and as we will see below in the study of the complete Hamil- 
tonian and its mean-field phase diagram the ^-fields will be 
the order parameter associated to the appearance of a topo- 
logically non-trivial phase. 

The contribution to the mean-field Hamiltonian H 2 mf 
arising from these off-diagonal terms, after transforming to 
Fourier space is given by 



^2MF diaS = ~ V 2 EE {2\Xa\M&a)n a - 3\ X a\ 2 ) 
k a 

with the momentum-dependent function 



(26) 



/ k (0) := cos(kai +0) +cos(ka 2 + 0) + cos(k(a! + a 2 ) -0). 

(27) 



Combining the two contributions 22 and (26) with the 
tight-binding term H t yields the expression for the com- 
plete mean-field Hamiltonian 

H 2 MF = E * f *J^2MF 00 * (k) 
k 

+ 3N A V 2 (\ X( p\ 2 + \x v \ 2 -(nl + nl)) (28) 



(m, n + 1) 



(m + l,n + 1) 



(m + l,n) 




FIG. 2. The rotational symmetry illustrated in this figure, in com- 
bination with translational symmetry, leads to the identification 
given by Eqs. (23) and (24). 



with 



^?2MF(k) 



V 2 (6n ( p-2\x ( p\f k {0<p)) 
-tAu 



V 2 {6n v -2\x v \f k i& v )))' 
(29) 



C. Complete Hamiltonian in the mean-field approximation 



By combining the results of the two previous subsections, 
we are in the position to formulate the mean-field approx- 
imation for the complete Hamiltonian (TJ. In momentum 
space it reads 



Hmf = X x F t (k)^f 2 MF(k) , i'(k) 

k 

+ 3N A (Vi {e - 7i 7ty) + V 2 [\ X ^\ 2 + \X V \ 2 ~ ml + n 2 ,))) 

(30) 



with the matrix elements of ^mf 00 given by 



(k) =3Vin v + &V 2 n,p - 2V 2 |^l/ k (0^ 
(k) =31^720 + 6V 2 n v - 2V 2 \x v \M® v ) 
Jtf™ (k) = [je™ (k)) * = - (f + f Vi) ai 



(31) 
(32) 
(33) 
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The Hamiltonian is readily diagonalized, yielding the zero- 
temperature free energy 

FIN A = 3V 1 $ 2 +3(Vi - 2V 2 )p 2 + 3V 2 {\x<p\ 2 + \% v \ 2 ) (34) 
-3{V l + 2V 2 )[n 2 -n) 

-v 2 ±f d 2 k{\ x <p\M©,p) + \x v \M®^) 

^ «/BZ 

^ [ d 2 k[{t+(v 1 f\A k \ 2 

Lr Jbz 



+ (3(Vi - 2V 2 )p + V 2 (|^|/ k (0 ) - \X V \M® V ))Y 



1/2 



where the expressions involving n a have been re-expressed 
in terms of the CDW order parameter p and the total 
fermion density n. Before starting with the calculation of 



the phase diagram, based on the free energy functional 34 , 
we add a few remarks: 

First, Eq. 34 compared to Eq. (18) shows that the repul- 
sive next-to-nearest neighbor interactions, which appear as 
the combination 3(Vi — 2V 2 )p in the integrand and in form 
of the constant 3(Vi — 2V 2 )p 2 in the free energy functional, 
effectively reduce the effect of the nearest-neighbor interac- 
tions, favoring that the system is in a phase with homoge- 
nous fermion density over the lattice, thereby tending to 
suppress the formation of a phase with CDW order. How- 
ever, for too strong next-to-nearest neighbor interactions 
V 2 , once (Vi -2V 2 ) < 0, one enters a regime in which in the 
employed two-lattice-site mean-field approximation the ef- 
fective interactions appear as if they were attractive®. In 
this limit the mean-field description used here obviously 
cannot be expected to reproduce reliable results. This de- 
ficiency can be dealt with by: (i) either working with a spa- 
tially enlarged MFT ansatz, which includes not a single two- 
site cell, but e.g. a six-site hexagon as the elementary spatial 
unit, as has been done in recent work by Weeks and Franz 62 . 
Such enlarged, more sophisticated ansatz correctly captures 
the repulsive character of the next-to-nearest neighbor in- 
teractions within the mean-field approximation, (ii) Alter- 
natively, one can stick to the spatially reduced two-site cell 
ansatz, and set for the described physical reasons the cor- 
rection due to V 2 in the 3(Vi~ 2V 2 )p and 3{Vi-2V 2 )p 2 terms 
in the free energy expression by hand to zero. This corre- 
sponds to not including the density-type decoupling chan- 
nel in the mean field approximation of the quartic V 2 term 
Hnnn (i-e. in suppressing the three terms resulting from the 
first line of Eq. (To) in the treatment of Hamiltonian (2l) that 
lead to the contribution (22)). As our main focus lies on the 
quantum simulation of the interacting fermion model on 
a cold atom platform, we follow the second approach (ii), 
which has also been adopted in previous work by Raghu et 

m 



Second, we remark that in the above derivation so far no 
assumptions have been made on the absolute values and 
phases of the order parameters Xcp and Xy ( see Section 



IIIB 



and Eq. (25)). To determine the phase diagram in the next 
section, we will adopt the restricted ansatz 



X<t> = +i\%\> X V = -i\xl 



(35) 



with fixed phases = ji/2, V - -n/2, and equally large 
absolute values. This par amet rization has been adopted 
previously in the literature ! 61 ! 62 ! Here we also choose this 
ansatz, also with the purpose to clarify a controversy about 
the numerical values of critical points in the phase diagram. 

Taking into account the described modifications, the free 
energy (34) reduces to 



1/2 



F/N A = 3Vif z + 3V lP 2 + 6V 2 \x\ 2 - 3(Vi + 2V 2 ){n 2 - n) 
g J^d 2 k [(f + ^i) 2 |4 k | 2 + (3 Vip - 2V 2 | X |g k (0)) 2 



(36) 



where the momentum-dependent function 
g k (0) := sin(kai+0) + sin(ka2 + 0) + sin(k(ai + a 2 )-0) (37) 
is related to / k (0) of Eq. (27) via g k (0) = / k (0-;z:/2). 

D. Mean-field phase diagram of the complete Hamiltonian 

Let us now determine the mean-field phase diagram. 



Minimization of the free energy 36 with respect to f , p and 



\X\ yields the coupled set of self- consistency equations 

* = 6l2 £ z d2fc (l^kl 2 (f+^l)^l)/W k (38) 



(39) 



P = ^of d 2 k{3V W -2V 2 \x\g k m)lW k 
2L Z Jbz 

III = ~-^^d 2 k [3V lP -2V 2 \ X \gk(0)) gk(0)/W k (40) 



with the dominator 



w k = (r + ^i) 2 l^ k l 2 + (3y 1 p-2y 2 lxlgk(0)) 2 



1 1/2 



(41) 



in the integrand of the three equations. Whereas in princi- 
ple all three equations have to be solved simultaneously, it 
is a good approximation to first treat the (relatively weak) 
effect of the renormalization of the bare hopping ampli- 
tude t, described by induced by the nearest-neighbor V\ - 
interactions independently, before dealing with the equa- 
tions for p and \%\, which are order parameters. 

In the region where interactions are weak and the system 
is in the SM phase (p = 0,x = 0), from Eq. 38 - or equiva- 
lently 19 - we find for the renormalization of the hopping 



amplitude, t' = t + ^Vi, 



6L 2 



jr. 

Jbz 



d 2 k\A k \ =0.262 Vi. 



(42) 



Next, we consider the limit with only repulsive nearest- 
neighbor interactions (Vy finite, V 2 = 0), for which the inter- 
actions favor the formation of a CDW phase. Indeed, from 
Eq. (20) we obtain 



7T=Tni d 2 k\AM\- x = 1.345 
Vic 2L Z Jbz 



(43) 
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FIG. 3. Phase diagram as determined from the free energy |36| 
obtained in MFT approximation. The blue lines correspond to 
second-order phase transitions from a semi-metallic (SM) phase 
at small interaction values to a CDW phase and a QAH phase. The 
red line describes a first order phase transition between the CDW 
and the QAH phase. The green color indicates the region of pa- 
rameter space, where a non-zero Chern number clearly indicates 
the non-trivial topological character of the QAH phase, as we will 
discuss in Sect.lIII Fl 



(44) 



and find that at a critical interaction strength of 
Vi c /t = 0.924 



in terms of the bare hopping amplitude, the system under- 
goes a second-order phase transition from the SM into a 
gapped CDW phase. Our results for the mean-field phase 
diagram, which we will derive and discuss step by step in the 
following, are summarized in Fig. [3] Figure[4]reveals details 
of the transitions between the different detected phases of 
this phase diagram, along three different paths in parame- 
ter space which are sketched in Fig.[4^. For instance, Fig.|4}D, 
which has been obtained by solving the self-consistency 



equation 19 for p for different V\ values, shows the buildup 
of CDW order p ^ for V\ > V\ c and makes manifest the 
second-order nature of the SM to CDW transition. 

Equivalently, we can consider the limit V\ = in which 
case we obtain 



Vic 



Y -r [ d 2 kg?{0)/\A k \= 0.845. 



3L 2 jbz 

For V\ = 0, where t' = t, the critical point 
V 2c /f = 1.183 



(45) 



(46) 



marks a second- order phase transition from the SM to a 
QAH phase, characterized by a non-zero order parameter 
\X\ ^ 0. This transition behavior is shown in Fig. [4]:, which 
has been obtained from the solution of Eq. (19) for different 
V2 values, and V\ - 0. 

Let us now comment on the earlier-mentioned contro- 
versy on the critical values of the second-order phase- 
transitions: our values for the renormalization of the hop- 
ping 1; and the critical interaction strengths V\ c and V20 at 



b. 

P 



® v 2 /t = c n 
1 




FIG. 4. Details of the continuous and first-order phase transitions 
between the SM, CDW and QAH phases found in the phase dia- 
gram of Fig.[3j. a) Different cuts in parameter space along the paths 
1 and 2 along the V\ and V-> axes (parts a) and c)) and along the 
path 3 from the CDW into the QAH phase (part d)). Green shading 
indicates parameter regions where a non-zero Chern number was 
found. 



which the transitions from the SM to the CDW and QAH 
phases set in, are in agreement with the values found by 
Weeks and Franz 62 . Our calculations confirm their critical 
values despite the fact that we work with a spatially reduced 
mean-field ansatz involving a simpler two-site basis cell as 
the elementary unit instead of a six-site basis cell 6 ^. Our 
two-site ansatz is, up to the effect of the renormalization 
of the hopping amplitude t' — t + £ V\ which we include in 
our treatment, equivalent to the ansatz chosen by Raghu et 
a 



(61 



If we neglect this hopping renormalization our Eq. 43 
yields t' IV\ C — ► tlV\ c = 1.345 from which we find a value 
Vic It - 0.743: this result is a factor of two smaller than the 
value found by Raghu et al, thus suggesting a problem of 
double counting of terms in the mean-field treatment as 
the origin of this discrepancy. We finally would like to re- 
mark that it is physically reasonable that the inclusion of the 
renormalization of the hopping due to the nearest-neighbor 
interactions V\ produces a larger value for the SM - CDW 
transition point {V\ c lt = 0.924 with hopping renormaliza- 
tion vs. Vic It - 0.743 without this renormalization): the 
Vi -induced increased hopping t' "acts in favor" of the SM 
phase and thereby pushes the transition point to the CDW 
phase towards larger Vi values. On the other hand, a moder- 
ate value of Vic 1 1 = 0.924 compared to the results found by 
Raghu et al. is good news from a quantum simulation per- 



spective (see Sect. IV below) as this value suggests that the 
CDW phase should be accessible already at more moderate 
values of effective nearest-neighbor interactions. 



Eq. 45 not only yields the critical value for V\ = 0, but 
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also describes the phase boundary between the SM region 
and the QAH phase for non-zero V\ -values, 



V 2c lt= 1.183 + 0.310 Vi/t. 



(47) 



Here, the non-zero slope of Vz c as a function of V\ is a re- 
sult of the renormalized hopping amplitude, which is in- 
creased for finite Vj values, thereby stabilizing the SM phase 
and pushing the QAH phase entrance point to larger Vz val- 
ues. On the other hand, for finite Vz # the phase bound- 
ary between the SM region and the CDW ordered phase re- 
sults to be a vertical line starting off at the critical value 
Vic/t = 0.9244 (see Eq. (44)). This behavior is rooted in the 



employed ansatz, where the 3(Vi - 2V2)p and 3(Vi - 2V->)p 



terms in the free energy are approximated b y 31/ 



IIIC 



p and 
Thus, 



3Vip 2 , respectively, as justified above in Sect 
within this ansatz any (weak) ^-dependence of the V\ 
driven SM to CDW transition is neglected. 




0.2 P 



0.1 P 



E. Free energy study 

Obtaining the phase diagram at arbitrary interaction val- 
ues, and in particular the study of the border between the 
CDW region and the QAH phase is more intricate as it re- 
quires the numerical solution of the coupled set of self- 
consistency equations (39) and 40 . From a technical point 



of view, numerically searching and finding solutions of the 
coupled self- consistency equations does generally not pre- 
clude the detected extremal points from corresponding to 
either (meta-stable) saddle-points or local minima in the 
free-energy landscape. It is therefore convenient to work di- 
rectly with the free energy functional: In contrast to the self- 
consistency equations based on local vanishing of a set of 
partial derivatives - the free energy captures the global fea- 
tures of the system. It allows to detect and study first order 
and second order phase transitions in a unified way. 

In particular, the presence of a first-order phase transi- 
tion between the CDW and the QAH phase as well as its 
underlying mechanism become especially transparent in a 
study of the free energy 36 . Figure [5] shows four contour 
plots which show the free energy landscape over the p-%- 
plane, for a fixed nearest-neighbor interaction strength of 
V\lt- 1.2 and increasing values of next-to-nearest neigh- 
bor interactions V2. For a small value V2I t - 0.1 the system 
is in the CDW phase, where the symmetry between <p- and 
i//-sublattice sites is spontaneously broken and the system 
resides in one of the two possible states with either posi- 
tive or negative relative fermion density imbalance +p on 
the two sublattices. These CDW order configurations cor- 
respond to the two global minima of the free energy, lo- 
cated on the p-axis (upper left part of Fig. [5j. For larger 
V2 values (upper right part of Fig. [5), a second pair of lo- 
cal minima of the free energy located on the j-axis appears. 
Despite the fact that the system is still in the CDW phase, 
the presence of these additional local minima at finite +\x\ 
values is a pre-cursor of the QAH phase and indicates that 
one approaches the transition into the topologically non- 
trivial phase. By further increasing V2 one reaches the criti- 



FIG. 5. Contour plots of the free energy (36) for a fixed value of 
V\lt= 1.2 and increasing values of V2. At small values of V2, the 
system is in the CDW phase. At increasing values for V2, a second 
pair of local free energy minima at finite \x\ values appears. At a 
certain V2 value these minima become global minima and the sys- 
tem undergoes a first-order phase transition from the CDW to the 
QAH phase. 



cal value V2 C (lower left part of Fig.[5j corresponding to the 
transition point at which the local minima with non-zero %- 
value become global minima of the free energy and the sys- 
tem enters the QAH phase (lower right part of Fig. [5j. This 
clearly indicates that the transition is of first order, which 
can also be seen from the abrupt jump of the p- and j-order 
parameter values across the transition, as shown in Fig.[4jd. 

This transition behavior with the associated change of lo- 
cal into global free energy minima and vice versa is also vis- 
ible in Fig.[6] which shows cuts of the free energy landscape 
along the p-axis (blue curve) and the |j|-axis (red curve). A 
comparison of the relative depth of the two types of minima 
allows one to precisely determine the numerical critical val- 
ues V20 which results in the red line denoting the first-order 
CDW to QAH transition shown in Fig. [3] One can show that 
for a fixed V\ interaction strength there exists only a single 
point V2c{V\) at which the depth of the two types of free en- 
ergy minima is the same, which in other words implies that 
the CDW-QAH transition is sharp and that a co-existence of 
the topologically trivial CDW phase with the QAH phase in a 
finite region of parameter space is precluded, at least when 
working with the mean-field free energy expression (36). 



F. Chern number characterization 

To unequivocally establish the topologically non-trivial 
character of the QAH phase, we have calculated the Chern 
numbers characterizing the topological nature of the sys- 



9 



tem in the various regions of the phase diagram, which for 
the Hamiltonian of interest in this work has - to the best of 
our knowledge - not been done so far. In contrast to the 
j-order parameter, which is a local order parameter, the 
Chern numbers are topological invariants and can be re- 
garded as non-local order parameters, which are only sen- 
sitive to the global properties of the system. They charac- 
terize the electronic states in the TMI phase of QAH type. 
As the ground state of the interacting QAH phase breaks 
time-reversal symmetry, this is described by a t/(l) Chern 
number^. 

Within the MFT approximation the Hamiltonian (30) is a 
quadratic form and can be readily diagonalized in momen- 
tum space, and its momentum-dependent energy eigenval- 
ues and eigenstates can be given analytically. We empha- 
size that despite this fact it is non-trivial to calculate the 
corresponding Chern numbers. The reason is that a nu- 
merical integration over the momenta of the Brillouin zone 
generically in one way or the other requires a discretiza- 
tion and the evaluation of the eigenfunctions and deriva- 
tives at a particular discrete set of momenta. It is a well- 
known problem by now 76 that if this discretization is done 
naively, the gauge invariance of the local curvature over the 
Brillouin zone is lost and the corresponding numerically ob- 
tained approximate Chern numbers fail to be an integer 
anymore, thereby loosing their intrinsic topological char- 
acter. Thus, to correctly obtain the Chern numbers which 
distinguish the topologically trivial and non-trivial phases 
of the system, we employ an efficient numerical algorithm 
which has been proposed by Fukui et at? 6 : The distinguish- 
ing feature of the method is that by construction of the al- 
gorithm the obtained Chern numbers are manifestly gauge- 
invariant and integer-valued, even when one works with a 
discretized Brillouin zone. The method only requires as an 
input the wave function of the system at such a finite set of 
points, without the need of working in a particular gauge or 
specifying gauge-fixing conditions. In addition, it produces 
the correct Chern numbers even for a fairly coarse-grained 
discretization of the Brillouin zone 76 . We refer the reader to 
Appendix[A]for more details on the numerical calculation of 
the Chern numbers in this work. 

Our numerical analysis confirms the topological QAH 
character of the % ji phase, as it unequivocally yields a 
Chern number value equal to one for the j-phase region 
and, as expected, vanishing Chern numbers for the topolog- 
ically trivial SM and CDW ordered phases (see Figs.|3]and|4]: 
& d) . We emphasize that while the order parameter % is local 
in the QAH phase, the Chern number is a non-local observ- 
able that provides additional information about the nature 
of this exotic phase. 



IV. QUANTUM SIMULATION WITH RYDBERG ATOMS IN 
OPTICAL LATTICES 

The dynamics of cold atoms in optical lattices is governed 
by Hubbard Hamiltonians 4 7 . The ability to control the un- 
derlying interactions and Hamiltonian parameters as well 




FIG. 6. Cuts of the free energy )36| as a function of p and \x\ for 
fixed V\lt and increasing values of V2 1 1. At small values of V->, the 
system is in the CDW phase. At larger values of V2, a second lo- 
cal free energy minimum appears. At a certain V2 value (around 
V%lt = 2.43) the latter becomes the global minimum and the sys- 
tem undergoes a first-order transition from the CDW to the QAH 
phase. 



as the dimensionality of the setup and the type of the lattice, 
in combination with sophisticated state preparation and 
measurement techniques, have turned cold atom experi- 
ments into a versatile platform for the study a wide range of 
complex many-body quantum systems, described by mod- 
els originating, e.g., from the field of condensed-matter 8 10 
or motivated from high-energy physics 89 . Specifically, as we 
will outline in detail in this section, cold fermionic atoms, 
loaded in an optical lattice and weakly laser- coupled to Ry- 
dberg states, provide a platform for a reliable quantum sim- 
ulation of the Hamiltonian of interest (see Eq. (TJ). 

Whereas the hopping dynamics in a two-dimensional 
honeycomb lattice geometry can be realized and con- 
trolled by a combination of sinusoidal optical trapping 
potentials^ 

the main challenge lies in the implemen- 
tation of the repulsive nearest- and next-to-nearest neigh- 
bor interactions of strength V\ and V-> between fermions 
residing at different lattice sites, at distances d nn = a and 
dram - V3a, respectively. Here we propose to use laser 
dressing of ground state atoms with strong and long-range 
Rydberg interactions to induce effective interatomic inter- 
actions over distances of several lattice sites. In particular 
we are interested in the parameter regime V\,V2 > t and 
V2 ~ V\, which is essential to be able to explore the pre- 
dicted quantum phases beyond the SM region of the phase 
diagram. 

Atoms laser-excited to high-lying Rydberg states^MJ in- 
teract strongly via long-range dipole-dipole or van-der- 
Waals interactions over distances of several fim. The in- 
teractions as well as the electronic lifetimes of atoms in 
Rydberg states are several orders of magnitude larger than 
for atoms in electronic ground states. One of the hall- 
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mark features of the strong Rydberg interactions is the 
"dipole blockade" effect, where large electronic level shifts 
inhibit the laser- excitation of more than a single atom to 
a Rydberg state within a certain blockade radius. In a 
many-body context, this mechanism gives rise to collec- 
tive effects and strongly correlated coherent many-body dy- 
namics in laser- driven atomic gased^Enil As suggested 
more than a de cade ag o 106 ! 107 * and recently demonstrated 
in the laborator y™ l 109 linternal state-dependent switchable 
Rydberg interactions can be exploited to realize neutral- 
atom entangling quantum gates, which are a key ingredi- 
ent for the development of scalable cold atom quantum in- 
formation processing architectures. In the field of quan- 
tum simulation the availability of (many-particle) Rydberg 
gates 110 111 allows one to realize a digital quantum simu- 
lation architecture, where coherent (as well as dissipative) 
time evolution of open many-body quantum syste ms is re- 
alized by stroboscopic sequences of gate operations 112 ^ 113 !. 
On the other hand, in the direction of analog quantum sim- 
ulation, the ability to tune the strength and form of the in- 
teractions by means of static and time-dependent fields al- 
lows one to engineer a broad range of Hamiltonians with ef- 
fective spin-spin interactions 1 14 " 1 18 i 



A. Rydberg dressing scheme and effective long-range 
interactions 

Following the analog quantum simulation approach, in 
the envisioned implementation, we aim at inducing effec- 
tive interatomic interactions by the application of "always- 
on" and global, i.e. spatially homogeneous CW-laser fields, 
without the requirement of single-site addressability, as re- 
cently achieved by several group a 119 l 120 [ To convey the main 
idea, we start by outlining first the ideal simulation scheme 
and comment afterwards on the range of validity, imperfec- 
tions and the experimental feasibility of our approach (see 
Sects. [IVBl and [r7cl below). 

We consider fermionic atoms loaded in a hexagonal opti- 
cal lattice, where - in addition to the external dynamics de- 
scribed by the hopping part of Hamiltonian {1} - the internal 
dynamics of the atoms is governed by a Hamiltonian 



H = Y 1 H i + Y 1 H iJ . 

i i<j 



(48) 



Here, in the co-rotating frame (and in rotating-wave approx- 
imation) the single-atom laser-driving term 



H i = (p.\r)(g\ i + h.c.)+A\r)(r\ i 



(49) 



accounts for coupling of atoms residing in a (meta) stable 
electronic ground state |g> to an excited Rydberg state |r> 
with a (two-photon) Rabi frequency O and detuning A, as 
depicted in Fig.[7^. Pairwise interactions between atoms ex- 
cited to the Rydberg state | r) are described by 



H l} = V i} \r){T\i9\r)(r\j. 



(50) 



We will focus on Rydberg s-states with repulsive van-der- 
Waals interactions Vij > 0, in which case the interactions 
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FIG. 7. a) Atoms residing in the stable electronic ground state |g> 
are coupled far off-resonantly (A » O) and in a red-detuned way 
(A > 0) to the meta-stable Rydberg state |r>. b) Schematics of the 
two-particle energy surfaces corresponding to the four two-atom 
electronic states; laser couplings between the electronic states are 
indicated by red arrows. The state | rr) with both atoms excited to 
the Rydberg state has a distance-dependent additional two-body 
energy shift due to repulsive van-der-Waals interactions Vn . 



Vij = c a \n 



~ a (with a — 6) are isotropic and fall off 



quickly with the interparticle distance. The pairwise inter- 
action potential curves as a function of the distance be- 
tween two atoms and the laser coupling between these 
curves are shown schematically in Fig.[7j3. 

For far off-resonant laser coupling, Q/A « 1, and the 
atoms initially residing in the electronic ground state |g), 
transitions from the ground state to states with one or more 
atoms excited to a Rydberg state are energetically not ac- 
cessible and strongly suppressed. Nonetheless, the laser 
coupling induces a weak admixture of Rydberg states to 
the ground state atoms (dressing), which thereby "inherit" 
some part of the Rydberg interaction character. To quanti- 
tatively describe the resulting effective interactions we de- 
rive the effective Born-Oppenheimer two-particle poten- 
tial surface for two Rydberg-dressed ground state atoms by 
means of standard (Van-Vleck-type) fourth-order perturba- 
tion theory 121 , whose validity is based on the small parame- 
ter fi/ A « 1. The resulting distance -dependent energy shift 
for two Rydberg-dressed ground state atoms is given by 



ft 4 

(A£)| gg> =2^3 



1 + 



2A 



(51) 



where we have subtracted the trivial, interaction- 
independent single-particle (fourth order in Q/A) AC- 
Stark-shift2(n 2 /A)(l-(ft/A) 2 )) of the two atoms in |g>. The 
effective two-atom interaction potential curve is shown in 
Fig.© 

In the weakly interacting limit (Vjj « A), at large inter- 
atomic distances, the energy shift |51[ reduces to 



CAE) 



\gg) ~ 



Vull 



2A 



(52) 



In this regime, both atoms can be partially excited to the Ry- 
dberg state, independently of the state of the other atom. 
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FIG. 8. Effective two-body interaction potential between Rydberg- 
dressed ground state atoms. The red curve shows the energy ob- 
tained from fourth- order theory (see Eq. |51|), whereas the blue 
curve is obtained from exact diagonalization of the two-particle 
Hamiltonian. At large distances, r > r c = (C6/(2A)) 1/6 (dotted, 
vertical line as a guide to the eye), Rydberg-dressed ground state 
atoms interact weakly according to 1/r 6 van-der-Waals interac- 
tions, which are reduced by a factor (H/A) 4 . In contrast, in the 
strongly interacting regime, at interatomic distances smaller than 
r c , the two-body interaction potential exhibits a plateau-like struc- 
ture, with a distance-independent energy shift 2 O 4 / A 3 . The black 
dashed lines correspond to the approximate analytical expressions 
for the two-atom potential in the weakly interacting and the strong 
blockade limit, as given by Eqs. |52) and |53[ , respectively. At even 
shorter distances 93 electronic wave functions of Rydberg atoms 
start to overlap and interactions are no longer described by van- 
der-Waals interactions. Parameters for this plots where chosen 

c 6 = a = ion. 



As a result, the interaction shift is in leading order given by 
the bare Rydberg interaction energy Vj-; = C a |r,- - tj \~ a , re- 
duced by the small prefactor (Q/A) 4 , which corresponds to 
the probability to find both dressed ground state atoms si- 
multaneously in the Rydberg state. On the other hand, in 
the strongly interacting limit (Vj ; « A) at small interparticle 
distances, where no more than one atom can be simultane- 
ously excited to the Rydberg, the energy shift becomes 



are naturally fixed by the lattice geometry, and the critical 
length scale r c , which in turn is determined by the choice 
of the electronic Rydberg state | r> and the tunable laser pa- 
rameters: To access the CDW phase it is sufficient to as- 
sure r c < d nnn , which implies V\ > Vz- Exploring the region 
where the existence of the QAH phase is expected, as well 
as its vicinity in parameter space, is more challenging: the 
crucial point is to realize a situation where pairs of nearest 
and next-to-nearest atoms interact strongly, i.e. are located 
within the blockade radius {d nn < r c and d nnn < r c ), where 
due to the plateau-like behavior of the potential curve (see 
Fig. [8j the interaction strengths are of comparable mag- 
nitude, V\ ~ Vt- Due to the continuous decrease of the 
Rydberg interactions as a function of the interparticle dis- 
tance inevitably some small, effective longer-range interac- 
tion terms (strength V3 for atoms separated by a distance 2a, 
V4 at a distance v7 a, etc.) will be induced. However, if these 
distances of at least 2a are larger than the critical radius r c 
the strength of residual longer-range interactions V3, V4, . .. 
fall into the weakly interacting regime, where they decrease 
rapidly with increasing distance (~ r -6 ) and are substan- 
tially smaller than the desired nearest- and next-to-nearest- 
neighbor interactions V\ and Vz. The question of the qual- 
itative and quantitative influence of these additional long- 
rang interaction terms on the phase diagram is interesting, 
and it has been partially addressed in recent work 62 , where 
the authors found that a third-nearest-neighbor repulsive 
interaction stabilized a QAH phase. However, a more de- 
tailed analysis of this question requires either a more so- 
phisticated (i.e. spatially extended) mean-field ansatz or ex- 
act diagonalization calculations, both of which lie beyond 
the scope of the present work. We note that inclusion 
of these terms might play a role in determining the pre- 
cise shape of the phase boundary, at which the system un- 
dergoes a transition into the non-trivial topological QAH 
phase. Whereas the mean-field analysis clearly indicates 
that strong next-to-nearest interactions are essential for the 
emergence of the QAH phase, the question whether Vj - Vz 
is sufficient to enter the non-trivial topological phase is a 
question to be ultimately answered by an experiment. 



O 4 ( A 



(53) 



Thus, up to a sub-leading correction, the energy shift is 
given by the distance-independent value of 2Q 4 /A 3 . The 
critical length scale r c at which the crossover from the 
plateau-like behavior with a universal value for the interac- 
tion shift at short distances to the weakly interacting regime 
takes place, is determined by the condition 2 A =; Vu, yield- 
ing 



r c = (C a /2A) 



I/O 



(54) 



The key idea to access the different regions of the phase 
diagram shown in Fig. [3] in the cold atom quantum simula- 
tion lies in controlling the relationship of the nearest- and 
next-to-nearest neighbor distances d nn and d nnn , which 



B. Validity of the effective two-body interaction potential 



In this subsection, we discuss the regime of validity for 



the effective two-body potential derived in Sect. IVA First of 



all, the off-resonant laser- coupling to the meta-stable Ryd- 
berg state I r> with a finite decay rate j r induces an effective 
decay rate je& f° r the dressed ground state atoms, which 
limits the description of the interactions in terms of the ef- 
fective two-body potential 51 to time scales t « y~J on 



which incoherent decay processes from the Rydberg state 
are negligible. In the limit j r « O, which is well-fulfilled for 
typical Rydberg states and moderately large Rabi frequen- 
cies Q, for the purpose of a reliable estimate one can simply 
include an imaginary part A|r)(r| — >■ (A - ij r l2)\r){r \ in the 
single-particle Hamiltonian Hi (see Eq. 49 ). For this non- 
hermitian Hamiltonian which includes the effect of decay 
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from the Rydberg state, the van Vleck-perturbation theory 
treatment then yields an effective decay rate per atom y e ff = 
(Q/A) 2 y r in leading order in (fi/A), both in the strongly and 
in the weakly interacting regime. At longer times such unde- 
sired decay processes from the Rydberg state set in, events 
which are typically associated to strong mechanical forces 
leading to atomic losses and heating of the dressed Rydberg 
gas, as recently studied in detail by Glatzle et alW^ 

In addition, it is a priori not clear that an effective descrip- 
tion in terms of effective pairwise two-body interactions is 
possible, and that three- and higher-order many-body in- 
teraction terms are negligible, in particular in the strongly 
interacting regime, where collective effects are expected to 
play an important role. As pointed out in previous woriP^H, 
the range of validity of the effective two-body description 
is determined by the interplay of the following conditions: 
on the one hand, due to the blockade radius r c the maximal 
density of Rydberg atoms that can be laser-excited in a two- 
dimensional lattice setup is given by n™^ = 1/r 2 . On the 
other hand, at large distances, the condition of weak dress- 
ing predicts a Rydberg atom density n^^ - (Q/A) 2 n, where 
n is the density of ground state atoms in the lattice. Thus, we 
find the condition nr% « (A/Q) 2 that has to be fulfilled. Us- 
ing the expression for the blockade radius (54) this amounts 
to the requirement 



n « 



2 A 



2-H 



= 23- 



vdW 1/3 



v 7/3 



o 2 c, 



1/3' 



(55) 



which can in principle always be fulfilled by choosing a suf- 
ficiently large detuning A at a given density n of ground 
state atoms. The ground state atom density is determined 
by the lattice geometry and the filling factor {n - 1/2 in 
the present case). Sixth-order in (Q/A) perturbation the- 
ory then shows that in the limit specified by 55 - and even 
in the strongly interacting regime, Vu » A - interactions 
are well-described by effective pairwise two-body terms ac- 



cording to Eq. 5 1 and that the strength of the leading-order 



corrections due to three-body interactions is of the order 
Q 6 /A 5 , i.e., at least two orders of magnitude smaller (in the 
parameter O/ A « 1) than the two-body terms. 



by means of the described Rydberg laser-dressing scheme, 
which are required for the implementation of the complete 
interacting model Hamiltonian (TJ. 

An optimal choice of Rydberg states and laser param- 
eters clearly depends on the specific experimental imple- 
mentation under consideration, the fermionic species used, 
and technical limitations such as available laser power and 
achievable two -photon Rabi frequencies for the coherent 
ground-to-Rydberg state laser excitation. However, let us 
provide a rough, though conservative estimate for a set 
of parameters, under which the requirements discussed in 
Sect. |IVA] are fulfilled and fundamental, non-technical lim- 
itations listed in Sect. IIVBI are small. We assume a real- 
space lattice spacing a of about 500 nrrP^l As discussed 
above, the exploration of the different phases of the model 
requires the critical radius r c (54) to be on the order of the 
lattice spacing a between neighboring sites. For bare van- 
der-Waals interactions of about V ~ 2n x 100 MHz for Ry- 
dberg s-states (angular momentum I - 0) with a principal 
quantum number around n — 30 at this distance^, this 
implies a detuning of A = V 12 — 2n x 50 MHz. Assuming 
a two-photon Rabi frequency of = 2n x 5 MHz yields the 
smallness parameter Q/A = 0.1 such the assumptions un- 
derlying the perturbative treatment are well-fulfilled. The 
resulting energy scale of the effective two-body interaction 
energy shift of dressed ground state atoms is (AE)\ gg ) ~ 
2fi 4 /A 3 = 2n x 10kHz. This energy corresponds to the V\ 
and V2 nearest- and next-to-nearest neighbor interaction 
strengths of Hamiltonian fil, thus being comparable to typ- 
ical tunnel rates t on the order of (a few) kHzP. 

Regarding imperfections: Rydberg s-states with a prin- 
cipal quantum number around n — 30 have typical radia- 
tive life-tim es Tr of a few tens of microseconds at room- 
temperatur c 193 ! 94 ! For x r ~ 30 /is we find an effective decay 
rate of y e ff = [Q.I A) 2 y r = 2n x 50 Hz, which is much smaller 
than the frequency scale corresponding to the effective in- 
teractions [AE)\gg). Finally, at half-filling of the lattice, for 
the present set of parameters the condition (55) is also ful- 
filled, such that interactions between laser-dressed ground 
state atoms are well-described by the effective two-body po- 
tentials, and three- and higher-order many-body effects are 
negligible. 



Experimental feasibility, Rydberg states and laser 
parameters for a quantum simulation 



V. CONCLUSIONS AND OUTLOOK 



Very recently, the simulation of (non-interacting) Dirac 
fermions with a quantum degenerate Fermi gas of 40 K 
atoms loaded into a tunable honeycomb lattice has been 
demonstrated experimentally 124 . On the other hand, sev- 
eral groups have recently achieved trapping and laser exci- 
tation of cold atoms in optical lattices into Rydberg states, 
and reported the observation of the Rydberg blockade in 
this setup^U^. If combined in a single experiment, these 
achievements directly enable the implementation of (i) the 
hopping dynamics of fermionic atoms in a two-dimensional 
hexagonal optical lattice potential, along with (ii) the engi- 
neering of the effective interactions over several lattice sites 



In this work, we have studied a system of interacting 
spinless fermions with nearest- and next-to-nearest neigh- 
bor repulsion on the honeycomb lattice and proposed a 
quantum simulation scheme of the model, based on cold 
fermionic Rydberg atoms in an optical lattice. In the first 
part of the manuscript, we have determined the phase di- 
agram of the model in a MFT treatment, yielding topologi- 
cally trivial SM and CDW ordered phases, as well as the exis- 
tence of a QAH phase for sufficiently strong next-to-nearest 
neighbor interactions. Beyond a quantitative comparison 
of our results with previous MFT studies, we have charac- 
terized the topological nature of the QAH phase by the nu- 
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merical calculation of a Chern number. In the second part 
of the work, we have proposed and worked out a scheme 
for an analog quantum simulation of this model, based on 
off-resonant laser-coupling of fermionic ground state atoms 
to electronically excited Rydberg states with strong van- 
der-Waals interactions. The proposed quantum simulation 
scheme is implementable with currently available experi- 
mental techniques, and it allows one in particular to access 
the regime of considerably large next-to-nearest neighbor 
interactions V2 — V\, which have been identified as one of 
the essential ingredients for the dynamical emergence of 
the QAH phase. 

The ideas that we have discussed here for Rydberg- 
dressed fermionic atoms can also be adapted to the 
platform of polar molecules, for which a rich tool- 
box for engineerin g of effective spin models has been 
developed 1 15 l 128 l 129 i and which have recently become avail- 
able in several laboratories^ 3 -^. 

For the quantum simulation of systems with topological 
properties with cold atoms or molecules in optical lattices, 
the use of engineered effective interactions to dynamically 
generate topologically non-trivial phases, as pursued in this 
work, constitutes a complementary route to the engineering 
of artificial gauge fields, which deserves further exploration. 
In this context, it will be important to understand the ef- 
fect of longer-range interaction terms, which are typically to 
some extent unavoidable in any physical implementation, 
on the stability of the encountered topological phases. On 
the other hand, the ample possibilies to engineer a rich va- 
riety of interaction potentials raises the question what other 
types of interaction- driven topological phases could be ex- 
plored in AMO quantum simulation setups, and how stable 
these fascinating novel quant um phas es are under realistic 
couplings to the environmen t 26 ! 131 -*^. 
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Appendix A: Numerical calculation of the Chern numbers 

Here we outline details of the numerical calculation of the 
Chern numbers based on the algorithm proposed by Fukui 
et a/f^-We apply the method to the mean-field Hamiltonian 



tonian reduces to 

H MF = E* t *)^2MF(k)#(k) 



with 



+ 3JV A (Vi(£ 2 -Vty) + 2 ^Il0l 2 ) (Al) 



3V l n v -2V 2 \x\Mnl2) -{t+V&Al \ 

-(t+V&Ak ■SV l n (l> -2V z \x\fv.(-nl2)} 



3V l n v + 2V 2 \ X \g k (.0) -{t+V^Al \ 
-(f+VtfMk 31^,-2^1x^(0)]' 



(A2) 



We work with a discretization of the Brillouin zone accord- 
ing to the grid k = [k x , k y ) with 

2n 

k x - l-Jivi, 

3 



2n 

iCy — — + fly S } 

3^ 

where the integer- valued indices n x , n y run as 

0< n x sC |47r/(3s)J, 
0< n y ^ L2?r/(v / 3s)J, 



(A3) 



(A4) 



such that the mesh of k- vectors covers the complete Bril- 
louin zone. The parameter s determines the degree of coars- 
ening of the mesh, and is fixed to s = 0.1 in our calculations. 
See also Appendix|B]for details on this and other equivalent 
parametrizations of the Brillouin zone. 

On this discrete set of points {k} we define the U(l) quan- 
tum link variables 



[/ M (k) = < M (k)|u(k+#>/iV 



(A5) 



Here, p, is a vector pointing to the next mesh point in the di- 
rection k x or k y , Nft a normalization constant of the hermi- 
tian product and | u(k)> is the normalized eigenvector corre- 
sponding to the lower-band eigenvalue £L (k) of matrix ]A2|, 



E-{k) = 3V 1 n-\{3V lP -2V 2 \x\gkmf+t' 2 \A k \ 2 ] , (A6) 

which can be readily obtained in analytical form. From the 
quantum link variables we then calculate the lattice field 
strength 

F xy Qt) =lnU x {k)Uy{k+ l k JU x (k+ 1* V^tfcOkT 1 . (A7) 



As emphasized in RefPH it is important and always possi- 
ble to choose a grid where for all k- vectors the eigenvectors 



in momentum space 30 , with the simplifying physical as- 



of Hamiltonian Al are well-defined. Inspection of Eq. A2 
shows that this can in the present case be guaranteed if a 
mesh is chosen where the points of the Dirac cones are ex- 
cluded (see Appendix|B]for explicit expressions for their lo- 
cations.) 

By evaluation of expressions JA5| and |A7| over the dis- 



sumptions discussed in Sect. |III C] under which the Hamil- crete set of points parametrizing the Brillouin zone, it is 
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FIG. 9. Left: Energy spectrum over the first Brillouin zone for van- 
ishing staggering potential ji = and t = 1. Right: Color-coded 
contour plot of the energy spectrum. 



then straightforward to compute the Chern number for the 
system at half-filling [n = 1/2), with the lower (upper) band 
completely filled (empty) , according to the sum over the dis- 
cretized Brillouin zone 



Ci : 



1 

2ni 



(A8) 



grid{ki 



Finally, we remark that the numerical algorithm is very 
efficient and sensitive in the sense that even for non-zero 
X values on the order of 10~ 10 , it yields a Chern number 
c\ — 1. As the precise value of % depends on the finite, nu- 
merical accuracy of determining the mean values for % from 



either minimizing the free energy 36 or solving the self- 
consistency equations (38) - 40 , we introduce a small, but 
finite cut-off below which we set non-zero ^-solutions to 
zero, to avoid spurious detection of non-zero Chern num- 
ber values, originating from finite numerical precision. 



Appendix B: Details on the Brillouin zone and domain of 
integration 



FIG. 10. The two blue- and gray-shaded triangles at the bottom 
of the hexagonal- shaped first Brillouin zone can be translated to 
the top of the hexagon to obtain an equivalent, rectangular- shaped 
area of integration. 




FIG. 11. Comparison of the Brillouin zone defined as the Wigner- 
Seitz unit cell of the reciprocal lattice (gray, hexagonal area) , and by 
an alternative, though equivalent definition as the area of the par- 
allelogram spanned by the vectors b\ and £?2 (blue-shaded area). 



In this appendix we provide explicit expressions and 
point out useful symmetry properties of various, equiva- 
lent parametrizations of the Brillouin zone of the hexago- 
nal lattice structure (cf. Fig.[l). In particular, we provide the 
parametrization which is used for the numerical calculation 
of the Chern number based on a discretized Brillouin zone. 

For the basis vectors of the real-space honeycomb lattice 
ai = (3/2, \/372) and a2 = (-3/2, v/3/2) (lattice spacing a set 
to one), the basis vectors of the reciprocal lattice are defined 
by a,.bj = 2^5;j and given by bi = (2nl3,2nl\/3) and b2 = 
(-27r/3,27z7\/3). The first Brillouin zone is a hexagon, for 
which the positions of the vertices are: 

4n In 2n 2n 2n 

'3\/3 ' 3 '3^3 ' 3 ' 3\/3 ' 
47T 2n 2n 2n 2n 
3^ ' 3 ' 3^3 ' 3 '3^3 ' 

The volume of the Wigner-Seitz cell of the real lattice, i.e. the 
area corresponding to one two-site basis cell, is 3\/3/2; thus 



the area of the domain of integration over the first Brillouin 
zone is L 2 = (2w) 2 /(3 v / 3/2) = 8^ 2 /(3\/3). Fig.[5]shows the 
two-band energy spectrum in the first Brillouin zone for 
vanishing staggering potential (5 = 0. 

Sometimes, in particular for comparison with results of 
related work in the literatu re 61 ! 62 ! it is convenient to move 
from the discrete sum over momenta to the integral over the 
Brillouin zone. Noting that 



£ 1 = N\, and / d 2 fc = L 2 , 



kEBZ J BZ 

this change can be realized according to 



E ■ 

kEBZ 



A/a 
L 2 



/■ 

Jbz 



d 2 fc. 



(B2) 



(B3) 



Furthermore, by making use of the translational in- 
variance of the Brillouin zone, one can change from the 
hexagonal-shaped Brillouin zone to a domain of integration 
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over a rectangular-shaped integration area, as indicated in 
Whereas the Brillouin zone, as shown in the left panel of 



Fig. 11 is defined in the standard way as the Wigner Seitz 



unit cell in the reciprocal lattice, it can alternatively also be 
defined as the parallelogram spanned by the vectors bi and 
b2 - as indicated by the blue-shaded area in Fig. 11 . In the 



present case of the hexagonal lattice the area is even more 
symmetric (a rhombus), since ||bil| = ||b2l|. 

Whereas the domains of integration of the rectangular- 
shaped and the parallelogram-like Brillouin zones are more 
suitable for numerical calculations, symmetry properties 
are most explicit in the standard definition (hexagonal- 
shaped Brillouin zone). 
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